1 Objective

The goal of this session is to get a clean macula densa dataset to work with.

2 Load Packages

options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320
## [1] 79456894976

3 Load Dataset

SO1 <- LoadSeuratRds(here("outputs", "so_merge_raw.rds"))

head(SO1@meta.data)

4 How many of these are actually MD cells

SO.markers <- FindAllMarkers(SO1, only.pos = TRUE)
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 5) %>%
    ungroup() -> top10
DoHeatmap(SO1, features = top10$gene) + NoLegend()

markers.to.plot1 <- c("Lrp2",         # PT
                      "Slc5a12",      # PT-S1
                      "Slc13a3",      # PT-S2
                      "Slc16a9",      # PT-S3
                      "Havcr1",       # Injured PT
                      "Epha7",        # dTL
                      "Cryab",        # dTL
                      "Cdh13",        # dTL1
                      "Slc14a2",      # dTL2
                      "Slc12a1",      # TAL
                      "Umod",         # TAL, DCT1
                      "Egf",          # TAL, DCT1,
                      "Cldn10",       # TAL
                      "Cldn16",       # TAL
                      "Nos1",         # MD
                      "Slc12a3",      # DCT
                      "Pvalb",        # DCT1
                      "Slc8a1",       # DCT2, CNT
                      "Aqp2",         # PC
                      "Slc4a1",       # IC-A
                      "Slc26a4",      # IC-B
                      "Nphs1",        # Podo
                      "Ncam1",        # PEC
                      "Flt1",         # Endo
                      "Emcn",         # Glom Endo
                      "Kdr",          # Capillary Endo
                      "Pdgfrb",       # Perivascular
                      "Pdgfra",       # Fib
                      "Piezo2",       # Mesangial
                      "Acta2",        # Mural
                      "Ptprc",        # Immune
                      "Cd74",         # Macrophage
                      "Skap1",        # B/T Cells 
                      "Upk1b",        # Uro
                      "Top2a"         # Proliferation
)
                      
DotPlot(SO1,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

DimPlot(SO1)

5 Subset out non MD cells.

SO2 <- subset(SO1, idents = c("4", "5", "6", "7"), invert = TRUE)

DimPlot(SO2)

SO2 <- SCTransform(SO2) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    FindClusters() %>%
    RunUMAP(dims = 1:30)
## Running SCTransform on assay: RNA
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 22835 by 11529
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 3 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 308 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 22835 genes
## Computing corrected count matrix for 22835 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 22.1136 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Klk1, Sostdc1, Wfdc2, Fabp3, Cox6c, Prdx5 
##     Ly6a, Foxq1, Wfdc15b, Atp1a1, Slc25a5, Krt7, Cox4i1, Cox5b, Ckb, Atp5g1 
##     Cox8a, Cldn19, Ndufa4, Atp1b1, Ggt1, Atp5b, Chchd10, Gm47708, Cox6b1, Uqcrq 
## Negative:  Pappa2, Malat1, Gm37376, Mir6236, Zfand5, Robo2, Neat1, Aard, Wwc2, Pde10a 
##     Itga4, Nadk2, Col4a4, Sdc4, Nos1, Col4a3, Sgms2, S100g, Ptgs2, Etnk1 
##     Ramp3, Irx1, Itprid2, Hnrnpa2b1, Srrm2, Tmem158, Ranbp3l, Bmp3, Camk2d, Gm24447 
## PC_ 2 
## Positive:  Ubb, Gm22133, Fth1, Mt1, Ldhb, Ftl1, Hspa1a, H3f3b, Fos, Eif1 
##     Junb, Cd63, Car15, Hspa1b, Prdx1, Rps24, Jun, Fxyd2, Mpc2, Cd9 
##     Actb, Mt2, Jund, Rpl7, Rps5, Mgst1, Btg2, Rps27a, Tmem213, Itm2b 
## Negative:  Malat1, Mir6236, Egf, Gm37376, CT010467.1, Umod, mt-Nd5, mt-Nd4l, Tmem52b, Nme7 
##     mt-Rnr1, Neat1, Slc12a1, Gm24447, mt-Nd6, mt-Nd2, mt-Co1, Atp1b1, Gm28438, Kcnq1ot1 
##     mt-Nd4, Lars2, Etnk1, mt-Rnr2, Wnk1, Dst, Foxq1, Sult1d1, Slc5a3, Atrx 
## PC_ 3 
## Positive:  Car15, Pappa2, Gm22133, Ldhb, Fth1, Mgst1, Itm2b, Cd63, Mpc2, Cd9 
##     Aard, Tmem213, Sfrp1, Ftl1, Fxyd2, Tmem59, Bsg, Ramp3, Prdx1, Clu 
##     Mcub, Tmem158, Irx1, Ppp1r1a, Sdc4, Ndfip1, Tmbim6, Tmem176a, Rpl9, Cldn10 
## Negative:  Fos, Junb, Jun, Egr1, Btg2, Hspa1b, Hspa1a, Atf3, Zfp36, Ier2 
##     Fosb, Klf6, Socs3, Klf2, Jund, Dusp1, Dnajb1, Rhob, Gadd45g, Malat1 
##     Gadd45b, Tsc22d1, Klf4, Ubc, H1f2, Cebpd, Ppp1r15a, Actb, Gm37376, Gm42793 
## PC_ 4 
## Positive:  Mt1, Gm8420, mt-Nd4l, mt-Co1, Gm8885, Rps21, Gm28438, mt-Rnr1, mt-Rnr2, Gm28437 
##     Apoe, Mt2, Rpl37, mt-Nd5, Rpl38, Gm28037, Gm11808, Gm3511, Rpl37a, Gm28661 
##     Rps28, Gm29216, Gm11353, mt-Cytb, Gm10925, mt-Nd2, Aebp1, Atp5md, Gm24447, Rpl39 
## Negative:  Pappa2, Umod, Egf, Aard, Tmem52b, Sult1d1, Tmem158, Car15, Mcub, Ptgs2 
##     Wwc2, Hsp90b1, Foxq1, Fth1, Dctd, Wnk1, Cd9, Gm22133, Iyd, Ptprz1 
##     Klk1, Sfrp1, Rpl15-ps2, Cdkn1c, Calr, Sec14l1, Atp1b1, Robo2, Clu, Tmsb4x 
## PC_ 5 
## Positive:  S100g, Gm8420, Tpt1, Rps21, Rps12, Pappa2, Rpl41, Rpl37, Gm8885, Rpl38 
##     Rpl37a, Mir6236, Rpl23, Rpl39, Gm11808, Gm6136, Actb, Rps27, Rpl23a, Gm3511 
##     Tmsb10, Rps26, Rpl30, Ppia, Tmem52b, Rpl5-ps1, Gm7206, Gm11353, Rpl34, Rps27a 
## Negative:  Hspa1a, Klk1, Hspa1b, Gm22133, mt-Rnr2, Car15, Atp1a1, Gm13340, Gm28437, Gm10925 
##     CT010467.1, mt-Cytb, Itm2b, mt-Nd4, Rpl15-ps2, Cd63, Clcnkb, Fth1, Gm28661, Ftl1 
##     Jun, Sfrp1, Gm29216, mt-Rnr1, Rpl13, Ptger3, Gpx6, Hspa8, Cldn10, Atp4a
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11529
## Number of edges: 377768
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7894
## Number of communities: 15
## Elapsed time: 1 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 08:58:06 UMAP embedding parameters a = 0.9922 b = 1.112
## 08:58:06 Read 11529 rows and found 30 numeric columns
## 08:58:06 Using Annoy for neighbor search, n_neighbors = 30
## 08:58:06 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 08:58:07 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuTO9ND/file150301d95057e
## 08:58:07 Searching Annoy index using 1 thread, search_k = 3000
## 08:58:08 Annoy recall = 100%
## 08:58:09 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 08:58:09 Initializing from normalized Laplacian + noise (using RSpectra)
## 08:58:09 Commencing optimization for 200 epochs, with 504322 positive edges
## 08:58:09 Using rng type: pcg
## 08:58:12 Optimization finished
DimPlot(SO2)

DotPlot(SO2,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead
## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead
## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b

6 Subset and Recluster

SO3 <- subset(SO2, idents = c("13"), invert = TRUE)

DimPlot(SO3)

ElbowPlot(SO3)

SO3 <- SCTransform(SO3) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:10) %>%
    FindClusters(resolution = 1) %>%
    RunUMAP(dims = 1:30)
## Running SCTransform on assay: RNA
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 22829 by 11472
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 278 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 22829 genes
## Computing corrected count matrix for 22829 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 19.84227 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Fabp3, Sostdc1, Foxq1, Klk1, Cox6c, Prdx5 
##     Wfdc2, Ly6a, Wfdc15b, Krt7, Cldn19, Ckb, Slc25a5, Atp1a1, Cox5b, Cox4i1 
##     Atp5g1, Cox8a, Ggt1, Ndufa4, Atp1b1, Gm47708, Atp5b, Gm44120, Chchd10, Gadd45g 
## Negative:  Pappa2, Malat1, Zfand5, Gm37376, Aard, Robo2, Mir6236, Wwc2, Pde10a, Neat1 
##     Itga4, Nadk2, Sdc4, Nos1, S100g, Ramp3, Col4a4, Irx1, Sgms2, Col4a3 
##     Ptgs2, Tmem158, Itprid2, Bmp3, Ranbp3l, Etnk1, Cdkn1c, Camk2d, Srrm2, Hnrnpa2b1 
## PC_ 2 
## Positive:  Fos, Junb, Jun, Hspa1a, Egr1, Hspa1b, Btg2, Atf3, Zfp36, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dusp1, Dnajb1, H3f3b, Tsc22d1, Rhob 
##     Gadd45g, Actb, Gadd45b, Mt1, Ubc, Ubb, Cebpd, Klf4, Mt2, H1f2 
## Negative:  Mir6236, Egf, CT010467.1, Gm37376, mt-Nd5, Malat1, mt-Nd4l, Slc12a1, Pappa2, mt-Rnr1 
##     mt-Co1, mt-Nd2, Nme7, Etnk1, Umod, mt-Nd4, Wnk1, Gm28438, Sfrp1, Gm24447 
##     mt-Nd6, Atp1b1, Robo2, Gm13340, Col4a4, mt-Nd1, mt-Rnr2, Kcnq1ot1, Sptbn1, Gm28661 
## PC_ 3 
## Positive:  Gm22133, Fth1, Ldhb, Ubb, Ftl1, Cd63, Car15, Mpc2, Cd9, Prdx1 
##     Mgst1, Eif1, Fxyd2, Mt1, Tmem213, Rps24, Itm2b, Bsg, Clu, Rpl7 
##     Rps5, Wfdc2, Rps27a, Mdh1, Rpl9, Rpl11, Aldoa, Ramp3, Rps7, Selenow 
## Negative:  Malat1, Mir6236, Gm37376, Egf, Fos, CT010467.1, Junb, Umod, Jun, Neat1 
##     Egr1, Nme7, Tmem52b, Gm24447, Atf3, Btg2, Zfp36, Fosb, mt-Nd5, mt-Rnr1 
##     mt-Nd6, mt-Nd4l, Klf6, Hspa1b, Kcnq1ot1, Socs3, Slc12a1, Ier2, Lars2, Klf2 
## PC_ 4 
## Positive:  Mt1, Gm8420, mt-Nd4l, mt-Co1, Gm8885, Gm28438, mt-Rnr1, mt-Rnr2, Apoe, Rps21 
##     Mt2, Gm28437, mt-Nd5, Rpl37, Aebp1, Gm28037, Rpl38, Gm28661, Gpx6, Gm10925 
##     Gm11808, Gm3511, Gm24447, Gm29216, Fgf9, Rpl37a, mt-Cytb, Egfl6, Atp5md, mt-Nd3 
## Negative:  Pappa2, Aard, Umod, Tmem52b, Egf, Sult1d1, Tmem158, Mcub, Foxq1, Ptgs2 
##     Wwc2, Car15, Hsp90b1, Dctd, Iyd, Ptprz1, Cd9, Wnk1, Fth1, Tmsb4x 
##     Cdkn1c, Ramp3, Gm22133, Rpl15-ps2, Gm44120, Calr, Clu, Sfrp1, Sec14l1, Col4a3 
## PC_ 5 
## Positive:  Hspa1a, Hspa1b, Gm22133, Klk1, Car15, Rpl15-ps2, Gm13340, mt-Rnr2, Atp1a1, CT010467.1 
##     Rpl13, Jun, Fth1, Itm2b, Ptger3, Cd63, Cldn10, Ftl1, Clcnkb, Hspa8 
##     Rpl10, Sfrp1, Gm10925, Mt2, Rpl9, Gpx6, mt-Nd4, Mt1, Atp4a, Lpl 
## Negative:  S100g, Gm8420, Gm8885, Rps21, Rpl37, Rpl41, Actb, Rpl38, Rpl37a, Rps12 
##     Tmsb10, Gm11808, Rpl39, Gm6136, Tpt1, Gm3511, Igfbp7, Tmem52b, Gm11353, Rps27 
##     Ppia, Rps26, Rpl5-ps1, Rpl23, Gm9843, Mir6236, Gm4332, Rpl23a, Pappa2, Gm12338
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11472
## Number of edges: 347067
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7654
## Number of communities: 14
## Elapsed time: 0 seconds
## 08:58:44 UMAP embedding parameters a = 0.9922 b = 1.112
## 08:58:44 Read 11472 rows and found 30 numeric columns
## 08:58:44 Using Annoy for neighbor search, n_neighbors = 30
## 08:58:44 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 08:58:44 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuTO9ND/file15030766e86ec
## 08:58:44 Searching Annoy index using 1 thread, search_k = 3000
## 08:58:46 Annoy recall = 100%
## 08:58:46 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 08:58:47 Initializing from normalized Laplacian + noise (using RSpectra)
## 08:58:47 Commencing optimization for 200 epochs, with 501470 positive edges
## 08:58:47 Using rng type: pcg
## 08:58:49 Optimization finished
DimPlot(SO3)

DotPlot(SO3,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead
## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead
## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b

FeaturePlot(SO3,"Pappa2")
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

FeaturePlot(SO3,"Nos1")

FeaturePlot(SO3,"Slc12a1")

6.1 Repeat Subsetting and reclustering

SO4 <- subset(SO3, idents = c("13"), invert = TRUE)

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:15) %>%
    FindClusters(resolution = 1) %>%
    RunUMAP(dims = 1:15)
## Running SCTransform on assay: RNA
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 22741 by 11431
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 319 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 22741 genes
## Computing corrected count matrix for 22741 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 21.52544 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Pappa2, Malat1, Gm37376, Zfand5, Mir6236, Aard, Robo2, Wwc2, Neat1, Pde10a 
##     Itga4, Nadk2, Sdc4, Nos1, Col4a4, S100g, Col4a3, Sgms2, Irx1, Ramp3 
##     Ptgs2, Etnk1, Tmem158, Itprid2, Srrm2, Hnrnpa2b1, Bmp3, Ranbp3l, Camk2d, Zbtb20 
## Negative:  Umod, Egf, Tmem52b, Sult1d1, Klk1, Sostdc1, Fabp3, Wfdc2, Cox6c, Prdx5 
##     Foxq1, Ly6a, Wfdc15b, Atp1a1, Krt7, Slc25a5, Cox5b, Cox4i1, Ckb, Atp5g1 
##     Cox8a, Cldn19, Ndufa4, Ggt1, Atp1b1, Atp5b, Gm47708, Chchd10, mt-Nd4l, Gm44120 
## PC_ 2 
## Positive:  Malat1, Mir6236, Egf, Gm37376, CT010467.1, Umod, mt-Nd5, Nme7, Tmem52b, Neat1 
##     mt-Nd4l, mt-Rnr1, Gm24447, Slc12a1, mt-Nd6, Kcnq1ot1, mt-Nd2, mt-Co1, Atp1b1, Gm28438 
##     Lars2, Etnk1, Wnk1, Dst, mt-Nd4, mt-Rnr2, Foxq1, Sult1d1, Slc5a3, Atrx 
## Negative:  Gm22133, Fth1, Ubb, Mt1, Ldhb, Ftl1, Eif1, Cd63, Car15, H3f3b 
##     Prdx1, Mpc2, Fxyd2, Cd9, Rps24, Mgst1, Tmem213, Hspa1a, Rpl7, Mt2 
##     Itm2b, Rps5, Bsg, Clu, Rps27a, Rps7, Rpl11, Aldoa, Rpl9, Mdh1 
## PC_ 3 
## Positive:  Pappa2, Car15, Itm2b, Mgst1, Cd9, Aard, Gm22133, Ldhb, Sfrp1, Mpc2 
##     Cd63, Tmem59, Fth1, Tmem158, Bsg, Mcub, Ramp3, Tmem213, Epcam, Prdx1 
##     Tmbim6, Ly6a, Gm13340, Clu, Hsp90b1, mt-Nd1, Rpl9, Tmem176a, Fxyd2, Ndfip1 
## Negative:  Fos, Junb, Jun, Hspa1b, Egr1, Btg2, Hspa1a, Zfp36, Atf3, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Gadd45g, Rhob, Tsc22d1 
##     Gadd45b, Ubc, Cebpd, H3f3b, H1f2, Klf4, H2bc4, Malat1, Actb, Gm42793 
## PC_ 4 
## Positive:  Pappa2, Umod, Aard, Egf, Tmem52b, Sult1d1, Tmem158, Car15, Mcub, Ptgs2 
##     Wwc2, Hsp90b1, Foxq1, Dctd, Fth1, Wnk1, Gm22133, Iyd, Cd9, Ptprz1 
##     Sfrp1, Rpl15-ps2, Cdkn1c, Klk1, Sec14l1, Calr, Clu, Ramp3, Robo2, Gm44120 
## Negative:  Mt1, Gm8420, mt-Nd4l, mt-Co1, Gm8885, Rps21, Gm28438, mt-Rnr1, mt-Rnr2, Gm28437 
##     Apoe, Mt2, Rpl37, mt-Nd5, Rpl38, Gm28037, Rpl37a, Gm11808, Gm3511, Gm28661 
##     Aebp1, Gm29216, Gm24447, Gm11353, Gm10925, mt-Cytb, Atp5md, mt-Nd2, Rpl39, Gm6136 
## PC_ 5 
## Positive:  Klk1, Hspa1a, Gm22133, Hspa1b, CT010467.1, mt-Rnr2, Gm13340, Rpl15-ps2, Car15, Atp1a1 
##     Rpl13, Itm2b, Gm10925, mt-Nd4, Fth1, Cd63, Gm28437, Ptger3, Ftl1, Rpl10 
##     Hspa8, Apoe, Clcnkb, Gpx6, Gm28661, Cldn10, mt-Cytb, Rpl9, mt-Rnr1, mt-Nd5 
## Negative:  S100g, Gm8420, Gm8885, Rps21, Pappa2, Rpl41, Rpl37, Rpl38, Actb, Rpl37a 
##     Tmem52b, Gm11808, Rps12, Tpt1, Gm6136, Rpl39, Gm3511, Ppia, Gm11353, Rps27 
##     Tmsb10, Aard, Gm9616, Gm7206, Rpl23, Gm9843, Rpl5-ps1, Rps26, Rpl23a, Gm12338
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11431
## Number of edges: 355957
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7629
## Number of communities: 16
## Elapsed time: 0 seconds
## 08:59:23 UMAP embedding parameters a = 0.9922 b = 1.112
## 08:59:23 Read 11431 rows and found 15 numeric columns
## 08:59:23 Using Annoy for neighbor search, n_neighbors = 30
## 08:59:23 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 08:59:23 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuTO9ND/file150307287c52f
## 08:59:23 Searching Annoy index using 1 thread, search_k = 3000
## 08:59:25 Annoy recall = 100%
## 08:59:26 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 08:59:26 Initializing from normalized Laplacian + noise (using RSpectra)
## 08:59:26 Commencing optimization for 200 epochs, with 475992 positive edges
## 08:59:26 Using rng type: pcg
## 08:59:29 Optimization finished
DimPlot(SO4)

DotPlot(SO4,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()
## Warning: Could not find Slc5a12 in the default search locations, found in 'RNA'
## assay instead
## Warning: Could not find Kdr in the default search locations, found in 'RNA'
## assay instead
## Warning: The following requested variables were not found: Slc13a3, Havcr1,
## Pdgfrb, Pdgfra, Piezo2, Upk1b

FeaturePlot(SO4,"Pappa2")

FeaturePlot(SO4,"Nos1")

FeaturePlot(SO4,"Slc12a1")

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:15) %>%
    FindClusters(resolution = 0.5) %>%
    RunUMAP(dims = 1:15)
## Running SCTransform on assay: RNA
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 22741 by 11431
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 319 outliers - those will be ignored in fitting/regularization step
## 
## Second step: Get residuals using fitted parameters for 22741 genes
## Computing corrected count matrix for 22741 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 22.5283 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Set default assay to SCT
## PC_ 1 
## Positive:  Pappa2, Malat1, Gm37376, Zfand5, Mir6236, Aard, Robo2, Wwc2, Neat1, Pde10a 
##     Itga4, Nadk2, Sdc4, Nos1, Col4a4, S100g, Col4a3, Sgms2, Irx1, Ramp3 
##     Ptgs2, Etnk1, Tmem158, Itprid2, Srrm2, Hnrnpa2b1, Bmp3, Ranbp3l, Camk2d, Zbtb20 
## Negative:  Umod, Egf, Tmem52b, Sult1d1, Klk1, Sostdc1, Fabp3, Wfdc2, Cox6c, Prdx5 
##     Foxq1, Ly6a, Wfdc15b, Atp1a1, Krt7, Slc25a5, Cox5b, Cox4i1, Ckb, Atp5g1 
##     Cox8a, Cldn19, Ndufa4, Ggt1, Atp1b1, Atp5b, Gm47708, Chchd10, mt-Nd4l, Gm44120 
## PC_ 2 
## Positive:  Malat1, Mir6236, Egf, Gm37376, CT010467.1, Umod, mt-Nd5, Nme7, Tmem52b, Neat1 
##     mt-Nd4l, mt-Rnr1, Gm24447, Slc12a1, mt-Nd6, Kcnq1ot1, mt-Nd2, mt-Co1, Atp1b1, Gm28438 
##     Lars2, Etnk1, Wnk1, Dst, mt-Nd4, mt-Rnr2, Foxq1, Sult1d1, Slc5a3, Atrx 
## Negative:  Gm22133, Fth1, Ubb, Mt1, Ldhb, Ftl1, Eif1, Cd63, Car15, H3f3b 
##     Prdx1, Mpc2, Fxyd2, Cd9, Rps24, Mgst1, Tmem213, Hspa1a, Rpl7, Mt2 
##     Itm2b, Rps5, Bsg, Clu, Rps27a, Rps7, Rpl11, Aldoa, Rpl9, Mdh1 
## PC_ 3 
## Positive:  Pappa2, Car15, Itm2b, Mgst1, Cd9, Aard, Gm22133, Ldhb, Sfrp1, Mpc2 
##     Cd63, Tmem59, Fth1, Tmem158, Bsg, Mcub, Ramp3, Tmem213, Epcam, Prdx1 
##     Tmbim6, Ly6a, Gm13340, Clu, Hsp90b1, mt-Nd1, Rpl9, Tmem176a, Fxyd2, Ndfip1 
## Negative:  Fos, Junb, Jun, Hspa1b, Egr1, Btg2, Hspa1a, Zfp36, Atf3, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Gadd45g, Rhob, Tsc22d1 
##     Gadd45b, Ubc, Cebpd, H3f3b, H1f2, Klf4, H2bc4, Malat1, Actb, Gm42793 
## PC_ 4 
## Positive:  Pappa2, Umod, Aard, Egf, Tmem52b, Sult1d1, Tmem158, Car15, Mcub, Ptgs2 
##     Wwc2, Hsp90b1, Foxq1, Dctd, Fth1, Wnk1, Gm22133, Iyd, Cd9, Ptprz1 
##     Sfrp1, Rpl15-ps2, Cdkn1c, Klk1, Sec14l1, Calr, Clu, Ramp3, Robo2, Gm44120 
## Negative:  Mt1, Gm8420, mt-Nd4l, mt-Co1, Gm8885, Rps21, Gm28438, mt-Rnr1, mt-Rnr2, Gm28437 
##     Apoe, Mt2, Rpl37, mt-Nd5, Rpl38, Gm28037, Rpl37a, Gm11808, Gm3511, Gm28661 
##     Aebp1, Gm29216, Gm24447, Gm11353, Gm10925, mt-Cytb, Atp5md, mt-Nd2, Rpl39, Gm6136 
## PC_ 5 
## Positive:  Klk1, Hspa1a, Gm22133, Hspa1b, CT010467.1, mt-Rnr2, Gm13340, Rpl15-ps2, Car15, Atp1a1 
##     Rpl13, Itm2b, Gm10925, mt-Nd4, Fth1, Cd63, Gm28437, Ptger3, Ftl1, Rpl10 
##     Hspa8, Apoe, Clcnkb, Gpx6, Gm28661, Cldn10, mt-Cytb, Rpl9, mt-Rnr1, mt-Nd5 
## Negative:  S100g, Gm8420, Gm8885, Rps21, Pappa2, Rpl41, Rpl37, Rpl38, Actb, Rpl37a 
##     Tmem52b, Gm11808, Rps12, Tpt1, Gm6136, Rpl39, Gm3511, Ppia, Gm11353, Rps27 
##     Tmsb10, Aard, Gm9616, Gm7206, Rpl23, Gm9843, Rpl5-ps1, Rps26, Rpl23a, Gm12338 
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11431
## Number of edges: 355957
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8189
## Number of communities: 8
## Elapsed time: 1 seconds
## 09:00:02 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:00:02 Read 11431 rows and found 15 numeric columns
## 09:00:02 Using Annoy for neighbor search, n_neighbors = 30
## 09:00:02 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:00:02 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuTO9ND/file150302fa2cb67
## 09:00:02 Searching Annoy index using 1 thread, search_k = 3000
## 09:00:04 Annoy recall = 100%
## 09:00:05 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:00:05 Initializing from normalized Laplacian + noise (using RSpectra)
## 09:00:05 Commencing optimization for 200 epochs, with 475992 positive edges
## 09:00:05 Using rng type: pcg
## 09:00:08 Optimization finished
SO.markers <- FindAllMarkers(SO4, only.pos = TRUE)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay:
## Gm13194, Gm10248, C430049B03Rik, C2cd4d, 1700022I11Rik, Bmp6, Gm6035, Gchfr,
## Gpc5, Gm3076, Gjb3, Gm13373, Gm7859, Il3ra, Gm3500

DimPlot(SO4)

7 Saving Seurat Object

save(SO4, file = here("jk_code", "JK_cleanMD.rds"))